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Abstract. We present results of numerical simulations of the propagation of 
ultra high energy cosmic rays (UHECRs) over cosmological distances, aimed at 
quantifying the statistical significance of the highest energy data on the spectrum 
and small scale anisotropies as detected by the AG AS A experiment. We assess the 
significance of the lack of a GZK feature and its compatibility with the reported 
small scale anisotropies. 

Assuming that UHECRs are protons from extragalactic sources, we find that 
the small scale anisotropies are incompatible with the reported spectrum at a 
probability level of 2 X 10~^. Our analysis of the AGASA results shows the power 
of combining spectrum and small scale anisotropy data in future high statistics 
experiments, such as Auger. 
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1. Introduction 

The problem of the origin of ultra-high energy cosmic rays (UHECRs) is a collection 
of several questions that will only be solved by high statistics observations around 
~ 10^" eV. The first question is whether there is a Greisen-Zatsepin-Kuzmin pP 
(GZK) feature in the spectrum of UHECRs due to photopion interactions with the 
cosmic microwave background (CMB). Photopion production is expected to generate 
a flux suppression around ~ 10^'^ eV possibly followed by a recovery. The statistical 
significance achieved by current experiments is not sufflcient to clarify whether the 
feature is present and what is its exact shape (see, e.g., but future data from the 
Pierre Auger Observatory |3] should resolve this issue. The second part of the problem 
is related to the identification of a class or several classes of sources that can accelerate 
particles to energies in excess of 10^" eV. The search for sources of UHECRs will likely 
be a long process requiring the identification of anisotropics in the arrival direction 
of the highest energy cosmic rays as well as multiparticle observations of high energy 
sources (see, e.g., Large scale anisotropics may indicate the generic host of 

sources while small scale anisotropics can eventually identify single sources directly. 
As small scale anisotropics start to become apparent with small number of events per 
source (e.g., doublets and triples), an estimate of the nmiibcr density of the sources 
becomes feasible (see, e.g. |S1C1)- For a clear source identification with single source 
spectral measuments, a future generation of experiments will be necessary. 

If astrophysical sources are responsible for the acceleration of UHECRs, the mean 
density of sources can be determined by analyzing small scale anisotropics (SSA) in the 
arrival directions that should become apparent as large data sets of the highest energy 
events become available. So far the only experiment to report SSA is AGASA [Hj and 
estimates for the number density of sources range from 10~^Mpc~'^ to ~ 10~'*Mpc~'^ 
with a large uncertainty |S] (see |7] for other estimates). In Ref. jS], a full numerical 
simulation of the propagation was performed accounting for the statistical errors in 
the energy determination and the AGASA acceptance with the conclusion that the 
AGASA small scale anisotropics is consistent with a density of sources ~ lO^^Mpc"'^. 
This source density corresponds to a luminosity per source of ^ lO'^^erg s^^ at energies 
E > IQ^^ eV and a best fit injection spectrum in the form of a power law has a spectral 
index 7 = 2.6 if the source luminosity has no significant evolution with redshift. For 
the case of strong source evolution, in [H] it was assumed that the luminosity varies as 
a function of redshift as L{z) cx (1 + z)* and the best fit slope at injection is 7 = 2.4. 
These fits are dominated by the energy regions with the largest number of events, i.e.. 
Sit E ^ 10^^ eV, such that a dominant contribution from galactic sources at these 
energies may significantly affect these results. The issue of the energy region where 
the transition from a galactic to extragalactic origin takes place is currently a subject 
of much debate. While the traditional view is that the ankle is indeed determined by 
such a transition, some recent work [HI llOj have proposed an alternative picture, in 
which the ankle is the natural result of the pair production process on cosmological 
propagation of extragalactic protons. The key issue in this class of models is the 
chemical composition, since a transition from heavy to light elements is predicted 
1^, 10_. A small contamination of the extragalactic component with iron nuclei would 
destroy the pair production dip and invalidate this model ^3 EI ■ If the ankle is due 
to this mechanism, the data at energies above and around the ankle are well fit by an 
injection spectrum in the form of a power law with slope 2.6 — 2.7 and a low energy 
cutoff that can be the result of the appearance of a magnetic horizon |11[I12|[T^ . 
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The results on small scale anisotropies, obtained for cosmic rays at energies above 
4 X 10^^ eV, are not significantly affected by a possible galactic contamination at lower 
energies. Moreover at such high energies the effect of the galactic magnetic field is 
expected to not affect the propagation in a significant way. 

The role of the intergalactic magnetic field is less constrained as different 
simulations give different estimates for the magnitude and spatial structure of these 
fields (see, e.g., ^|EI)- A homogeneous magnetic field spread through the whole 
universe would leave the spectrum of diffuse cosmic rays unchanged with respect to 
the case of the absence of magnetic field ^Sl- Here, we assume that intergalactic 
magnetic fields can be neglected for particles with energies above 4 x 10^^ eV. This 
assumption holds if magnetic fields in the intergalactic medium are less than ^0.1 nG 
if the reversal scale is 1 Mpc and the small scale anisotropies are evaluated on angles of 
~ 2 degrees . This magnitude field is compatible with observational bounds ^H] and 
detailed numerical simulations in JS] (however, see for different results). These 
magnetic field calculations are the result of the formation of the large scale structure of 
the universe. Other authors considered alternative possibilities, for instance, related 
to the early stages of evolution of galaxies ^HI , but in most cases the results are very 
model dependent and are affected by further evolution of magnetic fields during large 
scale structure formation. It is likely that UHECR observations will eventually allow 
us to measure the strength of the intergalactic magnetic field and get information 
about its structure (see, e.g., [T^On) V 

In the present paper, we address the significance of a detection of the GZK feature 
together with the identification of small scale anisotropies. We combine these two 
types of observations by AGASA in order to check their compatibility. In §2, we 
briefly describe the Monte Carlo code used to determine the spectra and the SSA. In 
§3, we discuss simulated spectra of AGASA and HiRes with regards to the presence or 
absence of a GZK feature. In §4, we calculate the two-point correlation function for 
different choices of the angular binning and combine the spectrum and SSA studies. 
We conclude in §5. 

2. The Monte-Carlo code 

The propagation of UHECRs is simulated using a Monte-Carlo code described in 
Plin]- We assume that UHECRs are protons injected with a power-law spectrum by 
extragalactic sources. The injection spectrum in taken to be of the form: 

F{E) cx E-^ cxp(-S/S„,ax), (1) 

where 7 is the spectral index and iSmax is the maximum injection energy at the source. 
We fix 7 = 2.6 and -Emax = 10^^'^ eV since these values reproduce well the experimental 
results in the lower energy, high statistics region ~ 10^^'^ eV to ~ 10^° eV Here, 
we focus on events above 4 x 10^^ eV, which are generated at z ^ 1, therefore, source 
evolution is only marginally relevant [S] . We simulate the propagation of protons from 
the source to the observer by including the photo-pion production, pair production, 
and adiabatic energy losses due to the expansion of the universe. The photo-pion 
production is treated as a discrete energy loss process and is simulated using the 
event generator SOPHIA |2T]. 

We use a constrained simulation where particles are propagated until the number 
of detected events above some energy threshold reproduces the experimental results. 
By normalizing the simulated flux by the number of events above a given energy, where 
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experiments have high statistics, we can estimate the fluctuations in the number of 
events above a higher energy where experimental results are sparse. In this way, we 
have a direct handle on the fluctuations that can be expected in the observed flux due 
to the stochastic nature of the photo-pion production and to cosmic variance. 

In choosing the initial particles to be propagated, we employ two different physical 
scenarios: a) a continuous distribution of sources; and b) a discrete distribution of 
point sources with a given density. For the continuous distribution a source distance 
is generated at random from a uniform distribution in a universe with JIa — 0.7 and 

= 0.3. In an Euclidean universe, every spherical shell contributes equally to the 
generated flux of particles. The flux from a source scales as r^^, where r is the distance 
between the source and the observer, while the number of sources between r and r + dr 
scale as r^, so that the probability that a given event has been generated by a source 
at distance r is independent from r. In a flat universe with a cosmological constant 
this is still true provided the distance r is taken to be 



where tg is the age of the universe when the event was generated, t^ is the present age 
of the universe, and R{t) is the scale factor of the universe. Once a particle has been 
generated, the code propagates it to the detector calculating energy losses and taking 
into account the energy and angular resolution of the detector. 

For a discrete distribution of point sources, we first generate the positions of 
sources with a given spatial density and then emit particles choosing randomly from 
the source positions. The source redshifts are generated with a probability distribution 
proportional to 

where r(z) is the r defined in Eq. |(5J and dt/Az gives the relation between time and 
redshift for a flat universe with a cosmological constant (see, e.g., the expression given 
in [2]). The source positions on the celestial sphere are chosen uniformly in right 
ascension and with a declination distribution proportional to cos(5.f 

Since we neglect the effects of the magnetic fields on the propagation, we ignore 
sources outside the experimental field of view and assign to visible sources a probability 
proportional to: 



where r(z)~^ takes into account the distance dependence of the solid angle and ijj{5) 
is the relative exposure of the experiment in the given direction. The positions 
for particle injection are chosen randomly from the assigned source positions and 
propagation proceeds as in the previous case. 

The last point worth stressing is the importance of the statistical errors in the 
energy determination. The statistical errors in the energy are accounted for in 
our simulation by assigning a detection energy chosen at random from a gaussian 
distribution centered at the arrival energy E and with width AE/E = 30%. As we 
show in the next section and in [J, the presence of a statistical error in the energy 
determination does affect the shape of the spectrum of UHECRs close to the GZK 
feature. This point was also made recently in j23j for the case of Log-normal errors in 
the energy determination. 

I The declination & is defined as being 0° on the equatorial plane and ±90° at the poles, so 
dn = cos SA&Aa. 




(2) 



r(z)-2^(<5). 



(4) 
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Figure 1. Effect on the spectrum of statistical errors on the energy 
determination. Outer plot: flux multiplied by the third power of the energy. 
Inner plot: flux. Dashed line: diffuse spectrum with no statistical error on the 
energy determination. Solid line: diffuse spectrum with 30% statistical error on 
the energy determination. 

3. Simulated spectra of AG AS A and HiRes 

There are several reasons to simulate the expected spectrum of UHECRs, rather than 
an analytical treatment, including the ability to introduce detector- related features 
which can affect the results. For instance, the energy dependence and zenith angle 
dependence of the acceptance of the experiment can be included as well as statistical 
and systematic errors associated with the energy detemination of each event. In order 
to illustrate the importance of this point, we calculated the spectrum of UHECRs 
from a uniform distribution of sources with and without statistical errors. In Fig. ^ 
we plot the expected fluxes (small square inside the main plot) with (solid line) and 
without (dashed line) 30% statistical error in the energy determination. The main 
plot illustrates the flux multiplied by the third power of the energy without statistical 
error (dashed line) , and with a 30% statistical error in the energy determination (solid 
line) . 

From Fig.^ it is clear that the effect of the statistical error on the flux multiplied 
by is to smear the shape of the GZK feature and to move its position toward higher 
energies. This effect is due to the steeply falling spectrum and to its slope changes, 
so that the probability of misplacing an event from a lower to a higher energy bin is 
larger than in the opposite direction. 

Results from the two largest exposure experiments, AGASA and HiRes, are 
claimed to be inconsistent. The fluxes from AGASA and HiResI are plotted in the 
left panel of Fig. |21 while the fluxes multiplied by are in the right panel. The flux 
figure suggests that despite the different techniques used by the two experiments, they 
appear to be in rather good agreement. There seems to be a systematic shift and some 
discrepancy at the highest energy bins. Their discrepancies are made more evident in 
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Figure 2. UHECRs spectra as measured by AGASA I2(il (circles) and HiResI 
1241 (squares). Left panel: normal spectra. Right panel: spectra multiplied by 
the third power of the energy. 




Figure 3. Histogram of the number of detected events with energy above 10^'' eV 
in 30000 realizations of the AGASA statistics. Dashed line: with no statistical 
error in the energy determination. Solid line: with 30% statistical error in the 
energy determination. The vertical line corresponds to the 11 AGASA events. 
The numbers above the last two bins are the number of events in those bins. The 
thin lines are gaussian fits to the histograms. 

the right panel, due to the amplification factor. At energies E > \Q™ eV, the HiRes 
spectrum shows a hint of a cutoff, while the AGASA spectrum continues unabated. In 
Ref. |5] we showed that a systematic overestimate of the AGASA energies by 15% and 
an underestimate of the HiRes energies by the same amount would in fact bring the 
two data sets in much better agreement in the region of energies below 10^" eV. At 
higher energies, the very low statistics of both experiments hinders any statistically 
significant claim for either detection or the lack of the GZK feature. 

To clarify the dangers of low statistics, we simulated the AGASA spectrum 
normalizing at 4x 10^^ eV (instead of 10^^ eV as in Ref. without the statistical error 
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on the energy determination and with a 30% error. We generated 30000 reahzations of 
spectra with 72 events above 4x10^^ eV hke the AGASA data and recorded the number 
of events with energies above 10^° eV in each reahzation. When the experimental error 
is taken into account, we obtained 3.5 ± 1.8 events above 10^° eV. In the case of no 
errors, the number is 2.5 ± 1.6. 

Assuming a gaussian distribution with the above numbers, we would incorrectly 
conclude that the 11 AGASA events above 10^° eV are a 4.2(t discrepancy, 
corresponding to a probability of 10"^. In Fig. 13 we plot the histogram of the 
number of events above 10^° eV in each realization. The dashed histogram is the 
result of the simulation without statistical errors, whereas the solid one is the result 
of adopting a 30% statistical error. The thin lines are gaussian fits to the histograms. 
As it is clear from Fig. 01 the distributions are not gaussian such that the naive 
estimates above of the discrepancy and of the probability are incorrect. Of the 
30000 realizations we simulated the number of spectra presenting 11 or more events 
with energy above 10^° eV is 18 for the simulations with statistical errors. This 
corresponds to a probability of 6 x 10~^ that in terms of gaussian errors would 
correspond to a discrepancy of about 3.2(t. A similar calculation was performed for the 
degenerate injection spectrum with 7 = 2.4 and redshift evolution of the luminosity 
L{z) oc (1 + zY- In this case the probability becomes 5.8 x 10"'*, basically identical to 
the previous case. 

In P], we also considered the effect of a possible systematic error in both the 
AGASA and HiRes experiments. To match the two spectra, we assumed that the 
AGASA energies had been overestimated by 15% and that the HiRes energies were 
underestimated by the same amount. In this case, the number of events above 10^" eV 
in the rescaled AGASA data set is between 6 and 7, while the number of events above 
]^q19.6 gY jg about 47. We simulated 30000 realizations of this statistics above 10^^'^ eV 
and counted the number of events in each realization with energies above 10^" eV. In 
this case, the number of realizations with 7 or more events above 10^° eV turned out 
to be 179, which corresponds to a probability of 6 x 10~^ or a discrepancy of 2.5cr 
in terms of gaussian errors. Since there are several events with energy very close to 
10^° eV we also considered the case of a number of events above 10^° eV equal to 6. 
Counting the simulated realizations with 6 or more events above 10^" eV, we obtained 
a probability of 0.024 or about 2a. 

These statistical considerations are the result of averages over large samples of 
realizations of source distributions, therefore one might wonder whether the spectra of 
UHECRs for those realizations that have a large number of events at ultra high energies 
are characterized by a pronounced GZK feature or are similar to the AGASA spectrum. 
To resolve this issue we plot, in Fig. ^ the spectra of some of the realizations that 
showed 11 or more events above 10^" eV. These spectra closely resemble the AGASA 
spectrum: all of them show no evidence of a GZK suppression, despite the fact that in 
the lower energy region they all fit the data quite well. This shows that an AGASA- 
like spectrum is not that improbable, even if the average cosmic ray spectrum can be 
expected to show a GZK feature. 

It is worth noticing that a few of the AGASA events are just a few percent above 
10^° eV, so that any small error in the energy assignment to one or more events would 
affect the statistical analyses discussed above. For instance, if just one of the events 
above 10^" eV had a real energy below 10^° eV, the probability of reproducing the 
AGASA result would become 2 x 10~^ (2.8cr). On the other hand, if one of the events 
below 10^° eV had actually an energy above 10^° eV then the probability would become 
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Figure 4. In these panels we plot 6 of the 18 simulations that have 11 or more 
events above lO^'' eV. The black crosses are the simulation results. The gray dots 
with error bars are the AGASA data superimposed for comparison. The error 
bars in the AGASA data are slightly shifted left to avoid covering of the black 



2 X 10^'' (3.6(t). So, changing the number of detected events above 10^° eV by just 
one unit changes the probabihties of reproducing the AGASA results by relatively 
large amounts. These results again indicate the excessively low statitics of events, 
and conclusions concerning the spectrum of UHECRs at these energies should not be 
considered as reliable. 

This argument applies as well to the case of the HiRes experiment. Although 
HiRes data appears to be consistent with the presence of the GZK feature, such 
a conclusion is also premature. To show this, we assumed that the actual flux of 
UHECRs is consistent with the spectrum measured by AGASA (which is usually 
considered to be inconsistent with the GZK feature) and selected events from such 
spectrum. Then each event is either accepted or rejected according to the HiRes 
aperture (the aperture of HiRes increases with energy while that of AGASA is 
constant at high energies). This procedure is repeated until the number of accepted 
events with energy above 4 x 10^® eV corresponds to the one measured by HiRes 
(namely 27). We then count the number of events with energies above 10^" eV in 
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each realization and calculate the probability of having one or less events (namely the 
HiRes number of detections) in that energy region. We find that such probability is 
~ 3%, corresponding to about 2cr in terms of gaussian errors. We also tried to vary the 
parameters related to the injection spectrum without appreciable changes in the final 
probabilities. Summarizing, HiRes is not unlikely to be observing events extracted at 
random from the AGASA spectrum. In other words, both AGASA and HiRes simply 
do not have the necessary statistics of events required to assess either the presence or 
absence of the GZK feature in the spectrum of UHECRs. 

4. Combining spectrum and small scale anisotropies 

If UHECRs are accelerated in astrophysical sources, then small scale anisotropies 
(SSA) should be observed at the highest energies. The level of SSA depends on 
the average density of sources and on the strength of possible magnetic fields in the 
intergalactic medium. Simulations of the formation of the large scale structure of the 
universe find that most of the volume of the universe is very weakly magnetized, and 
the typical deflections at energies above 4 x 10^^ eV are smaller than experimental 
accuracy [H) (however, see [TJ for a different conclusion). UHECR Astronomy should 
become possible at the highest energies starting with the detection of either small or 
large scale anisotropies. Observations of SSA can be used to infer the average density 
of sources and the typical luminosity of each source, an information which is not 
accessible through the spectrum of diffuse UHECRs alone. 

AGASA has reported the detection of SSA [5], whose significance has been 
questioned in [23|2H|- Even if not significant, the SSA data from AGASA can be 
used to exemplify the power of combining SSA and the spectrum. In the best 
fit density of sources based on AGASA SSA was estimated to be 10~^Mpc~^, 
although with large uncertainties (see also [?])■ In Fig. [SJ we plot the two point 
correlation function, defined as in [2] for different choices of the angular binning over 
which the multiplets are searched for (1 degree in the left panel, 2.5 degrees in the 
central panel, and 5 degrees in the right panel). The points with error bars are the 
expected (simulated) two point correlation function on the given angular scale for a 
homogeneous distribution of sources. It is clear that the statistical significance of the 
small scale anisotropies decreses significantly for angular binning on scales smaller or 
larger than ~ 2 degrees. 

In Section 3, we concluded that the measurements of the spectrum by AGASA 
and HiRes collaborations are not yet statistically significant around the GZK feature. 
In those calculations, we used the assumption of a uniform distribution of diffuse 
sources. We repeat these calculation for a discrete distribution of sources with density 
10-5 Mpc"^ With this source distribution, the average separation between sources 
is of the order of ^ 46 Mpc, such that the shape of the GZK feature should be more 
pronounced than for the uniform case. As a consequence, the disagreement of AGASA 
with the prediction of a GZK suppression should become more pronounced. In order 
to demonstrate this point, we simulate 50000 realizations of the AGASA statistics 
above 4 x 10^^ eV, in a scenario of discrete sources with density 10"^ Mpc"'^. In this 
case, we find 6 realizations out of the 50000 producing 11 or more events above 10^" eV 
and this corresponds to a probability of about 10"^ or 3.7cr. Clearly, the result is less 
(more) severe for source densities 10"'*Mpc~^ (lO"*' Mpc~^). 

This result is further illustrated in Fig. where we show a comparison between 
the spectra obtained for the AGASA statistics of events for a homogeneous distribution 
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Figure 5. Two point correlation function for AGASA data (liistograms) and 
Monte Carlo simulations (points with error bars). The latter are calculated for a 
homogeneous distribution of sources. 
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Figure 6. AGASA spectra multiplied by . The triangles are the result of our 
simulation in the scenario of continuous distribution of sources, the squares are 
the result of simulations with discrete sources with density of 10~^ Mpc~'^ and 
the dots are the AGASA data. 



of sources (black triangles with error bars) and point sources with density 10~^ Mpc~^ 
(blue squares with error bars). The circles are the AGASA data. 

In conclusion, the spectrum of AGASA is inconsistent (at the ^ 4(T level) with 
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Figure 7. In the above panels we plot the sky distribution of the simulated 
events in the 6 realizations that have 11 or more events above 10^" eV. 

a GZK feature generated by the source distribution implied by the reported SSA. 
These discrepancies could be more severe if the maximum energy at the accelerators 
is lower than that assumed above. In addition, the distribution in the sky of the 
11 (or more) events in these realizations are not necessarily similar to the AG AS A 
case. We show our simulated sets of events in Fig. [7| With the exception of the 
first map in the top left corner, all the realizations appear to have a higher number 
of multiplets when compared with AGASA observations. In other words, if UHECR 
sources are astrophysical proton accelerators, there is an inconsistency between the 
spectrum and small scale anisotropics of the AGASA experiment. If we require that 
both the spectrum and the sky distributions are correctly reproduced with a source 
density 10~^Mpc"^, the probabihty that this occurs in the AGASA data is ~ 2 x 10"^. 
This exercise shows that the combination of the spetral and anisotropy information 
is a powerful tool to estimate the density of sources of UHECRs and identify possible 
internal inconsistencies in the data. As the wealth of data expected from Auger 
becomes available in the next few years, this type of study should help the identification 
of UHECR sources. 
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5. Conclusions 

We investigated numerically the propagation of UHECRs using a Monte Carlo code 
originally developed in [2j and generalized to include astrophysical point sources in [HI . 
The code contains all the main particle physics aspects of the propagation, as well as 
several detector-related features (exposures, angle and energy dependent acceptances, 
statistical error in the energy determination, and limited statistics of events). We 
performed simulations constrained by the number of observed events above 4 x 10^^ 
eV and compared the number of events above 10^° eV observed by AGASA and HiRes. 
For AGASA we jointly analyzed the observed spectrum and small scale anisotropics. 
The combination of the two types of information strengthens the constraints on the 
presence of a the GZK feature or of clustering of events on small angular scales. 

In conclusion, we first presented the results of our simulations of the AGASA and 
HiRes statistics of events similar to the analysis in 5^ . We followed an alternative route 
and reached very similar conclusions. We simulated 30000 realizations with different 
scenarios and showed that the probability distribution of these events is non-gaussian. 
When expressed in terms of gaussian standard deviations the probability to record 
11 or more events with energies above 10^*^ eV corresponds to ^ 3cr. If a systematic 
energy decrease of 15% is introduced, the probability of a AGASA spectrum becomes 
~ 2%. We also found that the probability of finding a spectrum similar to HiRes 
(one or less than one event at energies above 10^" eV) is 3%. The two spectra cannot 
therefore be considered inconsistent with each other at a level better than 2a. 

We also analyzed the small scale anisotropics as detected by AGASA to stress the 
fact (already found in ^|) that the statistical significance of the clustering depends on 
the angular scale with which the data are binned. We clarified this issue by calculating 
the two point correlation function for the AGASA data and for the simulated data 
as obtained from a homogeneous distribution of sources. It is clear from this result 
that the angular binning chosen by AGASA is the one that optimizes the small angle 
clustering signal, but the statistial significance of such signal decreases for smaller or 
larger angular binnings. 

Finally, we combined the spectrum and the small scale clustering information. 
Choosing a source density of 10~^Mpc~'^ obtained by fitting the number of doublets 
and triplets of AGASA (assumed to be real), we recalculate the spectrum of UHECRs. 
As expected, for this source density the GZK feature is more pronounced than that 
obtained for a homogeneous distribution of sources and the discrepancy with respect 
to the AGASA spectrum is at the level of ~ Aa. Moreover, by simulating 50000 
realizations of the propagation, we obtained the arrival directions of the events in 
those realizations that show a number of events above 10^° eV equal to or larger than 
that of AGASA (11 events). Only 6 out the 50000 realizations fulfill the spectrum 
requirement and of the 6 realizations only one has a sky distribution of the events 
with resembles the AGASA one (no multiplets). We conclude that, if the sources of 
UHECRs are astrophysical proton accelerators, the reported small scale anisotropics 
are not compatible with the spectrum detected by the AGASA experiment at the level 
of 2 X 10~^ probability. In other words only 1 realization out of 50000 has a spectrum 
and a sky distribution of the events that are compatible with the AGASA data. 
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